Interfacial free energy of a hard-sphere fluid in contact with curved hard surfaces 
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Using molecular-dynamics simulation, we have calculated the interfacial free energy, 7, between 
a hard-sphere fluid and hard spherical and cylindrical colloidal particles, as functions of the particle 
radius R and the fluid packing fraction r\ — pa 3 /6, where p and a are the number density and 
hard-sphere diameter, respectively. These results verify that Hadwiger's theorem from integral 
geometry, which predicts that 7 for a fluid at a surface, with certain restrictions, should be a linear 
combination of the average mean and Gaussian surface curvatures, is valid within the precision of 
the calculation for spherical and cylindrical surfaces up to r\ ~ 0.42. In addition, earlier results 
for 7 for this system [Bryk, et al, Phys. Rev. E 68 031602 (2003)] using a geometrically-based 
classical Density Functional Theory are in excellent agreement with the current simulation results 
for packing fractions in the range where Hadwiger's theorem is valid. However, above rj » 0.42, 7(-R) 
shows significant deviations from the Hadwiger form indicating limitations to its use for high-density 
hard-sphere fluids. Using the results of this study together with Hadwiger's theorem allows one, in 
principle, to determine 7 for any sufficiently smooth surface immersed in a hard-sphere fluid. 

PACS numbers: 05.70.Np, 68.35.Md, 68.08.De 



The solid-liquid interfacial free energy, 7, is a central 
property governing a wide variety of technologically im- 
portant phenomena from crystal nucleation and growth 
to wetting. Because accurate and reliable experimental 
measurements of 7 are rare, much effort has been devoted 
in recent years to the development of atomistic simula- 
tion methods to determine this quantity for interfaces 
between coexisting solid and fluid phases [lrQ and for 
systems in which the solid is modeled by a static wall. 0- 
ll| These efforts have thus far been primarily restricted 
to planar interfaces; however, there are many physically 
relevant systems in which interfacial curvature is rele- 
vant, for example, in the formation of critical nuclei in 
nucleation [l2 - 14 or the solvation/wetting of hydropho- 
bic nanoscale particles. (l5- 17 [ There have been a num- 
ber of previous simulation studies that examine the effect 
of curvature in liquid- vapor interfaces HUH), but direct 
simulation studies on solid-liquid interfaces are lacking. 
In this work, we examine the dependence of 7 on the 
surface curvature for a hard-sphere fluid in contact with 
curved hard surfaces, specifically at spherical and cylin- 
drical colloidal particles. 

Konig, et aL[22( have recently shown that Hadwiger's 
theorem 23( from integral geometry puts severe restric- 
tions on the shape (curvature) dependence of the inter- 
facial free energy. In their analysis, the interfacial free 
energy of an object with a surface S is given by 

7 (5) = j + hH + KK (1) 

where h and k are constants depending upon the ther- 



modynamic state, but independent of the specific surface 
S. Here H and K are the averaged mean and Gaussian 
curvatures of S, defined as 
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where Ri(q) and R2(q) are the principal curvatures at 
each point q on the surface S, and A is the surface area. 
As discussed in Ref. [22j], the use of Hadwiger's theo- 
rem to determine the free energy of fluids at bounding 
surfaces, often referred to as "morphological thermody- 
namics" , is valid as long as the fluid/surface system satis- 
fies motion invariance, continuity and additivity. These 
conditions can break down for situations in which the 
bounding surface is small (i.e. on the order of the fluid 
particle size), for systems with long-range interactions, or 
for highly concave surfaces in which the fluid is confined 
in regions smaller than a few correlation lengths. For a 
fluid in contact with a spherical convex surface of radius 
R, H s = 1/R and K s = 1/R 2 , while for a fluid in con- 
tact with at convex cylinder of the same radius, one has 
H c = 1/2R and K c = 0. Much of the previous work on 
the curvature dependence focused on the first-order cur- 
vature correction, [24j referred to as the Tolman length 6 
which, for a convex spherical surface, would be given in 
terms of h by <$7o = —h/2. 



2 



The first theoretical treatment for the curvature de- 
pendence of 7 for a hard-sphere fluid at a hard wall was 
developed within the so-called Scaled Particle Theory 
(SPT),[lF " 
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|26| which is a theory of solvation that 
is based on an approximate determination of the work 
required to insert a spherical cavity (or hard-sphere so- 
lute) into a fluid. For this system, 7 scales trivially with 
T, that is, j(r),T) = 7* {rj)kT / a 2 , where 7* is the re- 
duced interfacial free energy, k is Boltzmann's constant, 
T is temperature and a is the hard-sphere diameter of 
the fluid. The SPT result for 7* for a hard-sphere fluid 
at a convex spherical surface of radius R is given by 



7spt(?7) 



3*7(2 + ??) 377 1 ln(l - rj) 1 

2tt(1 - 7]) 2 2tt(1 - 77) R 4tt R 2 
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where 77 is the packing fraction, defined in terms of the 
single particle density p as 77 — irpa 3 /6. (For simplicity 
in what follows we measure distance in units of a and 
energy in units of kT and drop the * superscript for re- 
duced units.) The SPT treatment is consistent with the 
Hadwigcr form of 7 (Eq. [TJ with h and k equal to the 
coefficients of 1 /R and 1/R 2 in Eq. @J respectively. Note 
that, the value of the interfacial free energy between a 
fluid and a wall depends upon the precise definition of 
the dividing surface that determines the volume of the 
bounding wall. In this work, we define the dividing sur- 
face to be coincident with the wall surface 

Using classical density functional theory (DFT) , Bryk, 
et al. 27[ have examined the scaling of the interfacial 
free energy and the excess interfacial adsorption for 
hard-sphere fluid at both spherical and cylindrical hard 
surfaces. This work, which utilized the Fundamental 
Measure Theory (FMT) version of DFT proposed by 
Rosenfeld[28|, determined 7 for this system as a func- 
tion of packing fraction and curvature with 1 / R ranging 
from to 0.5. The results for the spherical and cylin- 
drical geometries were fit to a polynomial of the form 
7(i?) = 7o + ai/_R + a 2 /-R 2 and were consistent with Had- 
wigcr's theorem in that the ratio of a\ for the spherical 
surface to that of the cylinder was 2. The predicted FMT 
values for h for both geometries are in good agreement 
with the SPT results (Eq.|4} at low packing fractions, but 
exhibit a negative deviation from the SPT curve of a few 
percent at the higher densities studied. More recently, 
Jin, et al.\ll\ examined the solvation free energy of var- 
ious ideal nanoscale non-spherical shapes (for example, 
cones, cylinders and prisms) using a hybrid Monte Carlo- 
DFT technique. They conclude that the use of "morpho- 
logical thermodynamics" based on Hadwiger's theorem 
to determine the solvation free energy (which includes 
both bulk and interfacial free energy components) for 



the non-spherical particles compares well with the DFT 
results. To date, there have been no determinations of 7 
from direct simulation for the hard-sphere fluid at curved 
surfaces with which to assess the accuracy of these DFT 
results. This is the principal goal of the present study. 

In this work, we make use of an adsorption equation 
derived using Cahn's extension (29| of the surface thermo- 
dynamics of Gibbs[3(|, namely, 
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where the excess interfacial volume per unit area, Vn, is 
defined by 
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where A is the interfacial V and N are the volume 
and number of particles, respectively, of a region contain- 
ing the interface and Vf and Nf are the correspon ding 
quantities for a region entirely within the bulk fluid. 
This adsorption equation can be shown to be equivalent 
through a change of variables and Maxwell relations to 
the usual Gibbs adsorption equation 



dp 
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where r^r is the excess interfacial number of particles per 
unit area at the interface. Eq. [5] can be integrated with 
respect to pressure to give 



j(P) = 7P=o+ / v N {P)dP 



(7) 



where P is the pressure. For the hard-sphcrc/hard-wall 
case, 7p=o = 0. Eq. [7] was recently used to determine 
7 for a hard-sphere fluid at a flat wall[ll|. The excess 
interfacial volume can be related to the density profile 
p(r) by 
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where i = 1 and 2 for the cylindrical and spherical ge- 
ometries, respectively. 

To determine 7 for this system using Eq. [71 we use 
molecular-dynamics (MD) simulation to calculate vn as 
a function of bulk packing fraction r\ = pa 3 /6 for a hard- 
sphere fluid in contact with spherical and cylindrical col- 
loidal particles of radius R varying from 10 down to 0.5, 
over a range of r\ from 0.03 to 0.49, the upper limit be- 
ing the fluid packing fraction at freezing. For the sim- 
ulations, we use the algorithm of Rapaport [3~fl | and use 
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Eq. [5] directly to determine vn in the simulations. For 
additional simulation details and plots of v(rj;R) see the 
Supplemental Information. [32j 

Because we find as a function of the packing frac- 
tion, not the pressure, we transform Eq. [7] to give 
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(9) 



To obtain the derivative of P with respect to rj we use 
the KLM-low equation of state (EOS). This EOS has 
been shown to give five decimal place accuracy in the 
pressure even at hi gh density when compared to high- 
quality simulations, [34| so any errors introduced by its 
use are much smaller than the statistical error in the sim- 
ulation data. To reduce the numerical integration error 
in evaluating Eq. [9] we subtract from the integrand the 
corresponding value obtained from the SPT (Eq. [4} - the 
corresponding SPT expression for the cylindrical surface 
is obtained from Eq. @] by setting the coefficient of the 
1/R 2 term to zero and dividing the coefficient of the 1/R 
term by two. Accurate calculation of the excess interfa- 
cial volume v at very low densities is difficult because of 
sampling issues; however, this quantity can be calculated 
exactly in the limit p — > 0. In this limit, we have 
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where a\ and <Z2 are constants equal to 1/4 (1/8) and 
1/24 (0), respectively, for the spherical (cylindrical) wall. 

For small values of R, the statistical errors in vn are 
larger for similar simulation lengths because the number 
of particles near the wall is relatively small. For the spe- 
cial case of R = 0.5 where the spherical surface is iden- 
tical in size to the fluid particles, high precision results 
can be generated by recognizing that the density profile 
p(r) is equivalent to pfg(r), where g(r) is the radial dis- 
tribution in the bulk hard sphere fluid. For this value 
of R, wc use an MD simulation for a bulk hard-sphere 
fluid to calculate g(r) from which vn can be determined 
with high precision. Alternatively, by replacing p(r) with 
Pfg(r) for R = 0.5, we have from Eq. |8]that 
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v N (R = 0.5) = 4 / r 2 [l - g{r)]dr 

J 1/2 



(11) 



4irr 2 [g(r) - l]dr - 1/6 



Using the compressibility equation from liquid state 
physics, [35| the integral in the previous equation can 
be replaced with a term dependent upon the isothermal 



compressibility kt = p 1 (dp/dP)T to yield 

v N (p; R s = 1/2) = --{kTn T - p- 1 ) - 1/6 (12) 

TT 

Using an equation of state, the isothermal compressibil- 
ity Kt can be determined analytically as a function of rj 
allowing for an analytical calculation of vn, and thus 7 
for R=l/2. 

The calculated values of 7 from our simulations for the 
spherical and cylindrical walls are plotted as functions of 
r\ in Fig. [TJ for values of R ranging from the planar wall 
[R = 00) to 0.5. At all densities, these figures show that 
at fixed 77, the intcrfacial free energy 7 is a monotonically 
increasing function of 1/R. 
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FIG. 1. Color online. Top panel: Interfacial free energy be- 
tween a hard-sphere fluid and a hard cylindrical colloidal par- 
ticle as a function of rj for several values of the radius, R. 
Bottom panel: Same as the top panel except for the spherical 
surface. The inset shows the value of the excess volume, vm, 
for the spherical surface as a function of packing fraction for 
R = 00, 5 and 0.5. The symbols in this panel and in the inset 
are as indicated in the top panel legend. 

To test the validity of Hadwiger's theorem for this sys- 
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tern, we fit the data shown in Fig.[T]to Eq.[T]to determine 
the constants h and k for each 77. The fit was performed 
using standard weighted quadratic (for the spherical case 
where H = 1/R and K = l/R 2 ) and linear (for the cylin- 
drical case where H = 1/2R and K = ) least-squares 
regression. The weights in the regression were equal to 
the inverse of the statistical variance of the data points. 
The results of these fits are shown in Fig. [2] for packing 
fractions up to 0.42. A table of the fitted values for 70, 
h and k is included in the Supplemental Information. |32j 
Also, shown in Fig. [2] are the results for h and k from the 
SPT expression Hand from the Rosenfeld DFT.[27| The 
DFT results agree remarkably well with the simulation 
results for the values of 77 shown in Fig. [2] Most notable 
in these results is the fact that the estimates for h for the 
spherical and cylindrical geometries shown in Fig. [5] are 
in excellent agreement, as predicted by Hadwiger's the- 
orem. The SPT prediction agrees well with the simula- 
tion results at low packing fractions, but underestimates 
h and slightly overestimates n at packing fractions above 
about 0.25. Note that, as mentioned earlier, the value of 
7 and, by extension, the Hadwiger coefficients 70, h and 
k will depend upon the choice of dividing surface; how- 
ever, the Hadwiger form (Eq. [TJ is remains valid, albeit 
with modified coefficients. }l6| 

Above 7/ = 0.42, however, significant deviations from 
the Hadwiger form were found. To quantify this devia- 
tion, we performed a cubic and quadratic weighted least- 
square regression with respect to 1 / R on the sphere and 
cylinder data, respectively. The results of these calcu- 
lation are plotted in Fig. [3] At packing fractions below 
0.42, the cubic coefficient (03) for the spherical surface 
is zero within the estimated statistical error, consistent 
with Hadwiger's theorem, but this coefficient diverges 
quickly from zero at packing fractions above 0.42, in- 
dicating a significant breakdown of the Hadwiger form 
at packing fractions approaching the freezing density 
(77 = 0.492). The situation is similar for the quadratic 
coefficient (02) for the cylindrical case (Fig.[3]top panel)- 
there is one value of a-i at 77 0.25 that is marginally dif- 
ferent from zero outside the error bars, but given its 
marginality and the small value of the coefficient, we do 
not view this deviation as significant. The origin of the 
divergence of the non-Hadwiger coefficients at high pack- 
ing fraction can be seen in the inset in Fig. [1] in which the 
excess volume, uat, is plotted for the spherical surface as 
a function of 77 for R/a = 00, 5, 1 and 0.5. At high 77, the 
excess volume v n for the planar wall (R = 00) exhibits 
a downward curvature that is not present in the other 
values of R, making it impossible to fit vn(R) purely 
as a quadratic polynomial in 1/R. It is possible that 
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FIG. 2. Color online. The calculated values of h and k for the 
spherical and cylindrical walls. Filled symbols are the results 
from the present MD simulations and the open symbols are 
the values calculated from the White-Bear DFT. The solid 
and dashed lines are the values of h and n from Scaled Particle 
Theory (SPT). Note that, the symbols for the various values 
of h are difficult to resolve from one another because they are 
nearly coincident. 



this anomalous decrease in Vn for the planar wall is con- 
nected to the prefreezing transition that has been ob- 
served for the hard-sphere fluid at a hard planar wall 
at high packing fraction (Tol 36, 37 1, and further analysis 
is underway. These results also indicate that the errors 
in the SPT at the intermediate packing fractions are not 
likely to be due to the omission of higher-order terms (be- 
yond the quadratic for the sphere) in 1/R, as has been 
proposed (see Ref. [38| for a review). The results here 
put severe upper bounds on the magnitude of any such 
higher-order terms. As an example, our calculated value 
for 7 at 77 = 0.25265 and R = 1.0 is 0.6595(6). The cor- 
responding SPT value from Eq. [His 0.6711 - a difference 
of 0.0116(6), which is more than an order of magnitude 
larger than the estimated upper bound to the contribu- 
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linear combination (h and k) as functions 77. For these 
packing fractions, the Hadwiger form for 7 (Eq. [1} is 
shown to be valid for both convex spherical and cylin- 
drical surfaces even down to R = 0.5a, where the radius 
of the bounding surface is equal to that of the fluid par- 
ticles. In addition, earlier DFT results based on Rosen- 
feld's fundamental measure theory are found to be in 
excellent agreement with the results of the current simu- 
lations. At higher packing fractions (77 > 0.42), however, 
significant deviations from the Hadwiger form are ob- 
served. The results of these simulations can serve as a 
useful reference model to determine the thermodynamics 
of the hydrophobic solvation of nanoscale particles. 
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FIG. 3. Bottom panel: Plot of the first non-Hadwiger co- 
efficient (cubic, az) of 1/7? in the curvature expansion of 7 
for the hard-sphere fluid at a spherical colloidal particle as 
a function of packing fraction 77. For clarity, the inset shows 
the data for 77 < 0.42 on a smaller scale. Top panel: Same as 
bottom panel, but for the cylindrical surface. 



tion of the cubic term CI3 at this packing fraction from 
Fig. [3] Note that largest contribution to the difference 
comes from 70, the value for the planar wall. 

In summary, we have calculated the curvature depen- 
dence of the interfacial free energy between hard-sphere 
fluid and hard spherical and cylindrical colloidal parti- 
cles of varying radius, R. These calculations have im- 
portant applications in determining the solvation free 
energy of nanoscale objects. From our simulation re- 
sults, we have verified that the predictions of Hadwiger's 
theorem, [22|, 23 1 which predicts that 7 for fluids at a 
bounding surface, under certain restrictions, is com- 



pletely described by a linear combination of the aver- 
age mean and Gaussian curvatures, arc valid for the sys- 
tems studied at low to moderate fluid packing fractions 
(77 < 0.42) and have calculated the coefficients of this 
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